Efficient algorithm for multi-qudit twirling for ensemble quantum computation 
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We present an efficient algorithm for twirling a multi-qudit quantum state. The algorithm can 
be used for approximating the twirling operation in an ensemble of physical systems in which the 
systems cannot be individually accessed. It can also be used for computing the twirled density 
matrix on a classical computer. The method is based on a simple non-unitary operation involving 
a random unitary. When applying this basic building block iteratively, the mean squared error of 
the approximation decays exponentially. In contrast, when averaging over random unitary matrices 
the error decreases only algebraically. We present evidence that the unitaries in our algorithm can 
come from a very imperfect random source or can even be chosen deterministically from a set of 
cyclically alternating matrices. Based on these ideas we present a quantum circuit realizing twirling 
efficiently. 

PACS numbers: 03.67.-a,03.67.Lx,02.70.-c 



I. INTRODUCTION 

Twirling was first introduced for bipartite systems in 
Refs. [I], Q in the context of entanglement purification 
and still appears as part of various quantum information 
processing protocols 0,0, IE]- For example, in the single 
party case, twirling makes possible to obtain the average 
gate fidelity of a positive map [6|. Later twirling was 
generalized to multipartite systems: For a given density 
matrix p the twirled state is defined as 0, d, Q 



Pp:= 



ueu(d) 



u® N P (u® N )Uu, 



(1) 



where U(d) is the group of d-dimensional unitary matri- 
ces, iV is the number of qudits, and dU is the normalized 
Haar measure over U{d). For the bipartite case one can 
also consider twirling defined as 



U <8> U* p(U ® U*)UU, 



(2) 



ueu(d) 



where '*' denotes element-wise complex conjugation. 
States obtained from P and Pi so are called Werner states 
and isotropic states, respectively. Isotropic states are 
quite useful in quantum information processing: They 
are the maximally entangled state mixed with white 
noise. While in this work we focus on P, our results 
generalize trivially to the computation of P; so . 

The importance of twirling in the multipartite case 
is that it transforms a general mixed state into a state 
that can be characterized with only a few parameters 
0> E3, El ■ Entanglement of formation is known for bi- 
partite isotropic states [H, [H[ and necessary and suffi- 
cient conditions for the entanglement of tripartite Werner 



states are also known ||. Therefore, since twirling can- 
not increase any entanglement monotone, if we can ex- 
perimentally twirl a state, we can simplify the estimation 
of its entanglement properties. 

Moreover, twirling also appears in various calculations 
in quantum information science (e.g., it is used to define 
a family of quantum states in Ref. Integrals over 

U{d), similar to twirling, appear in many areas of physics 
[15j | . In particular, the computation of integrals of the 
form 



Ui 1 j 1 Ui. 



■U, 



XU kl hU k2l2 ...U knln )*dU, (3) 



ueu(d) 
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are needed. Such integrals for the m = n case can 
straightforwardly be obtained from twirling appropri- 
ately chosen density matrices. Twirling is also closely 
related to unitary t-designs which have raised interest 
recently [II Hi- 
lt seems straightforward to implement twirling: One 
has to apply a random multilateral unitary rotation to 
each copy of a state, and then average over the ensemble. 
However, problems quickly arise when considering prac- 
tical implementations. Applying different random rota- 
tions to different systems of the ensemble requires that we 
are able to access the systems individually. In practice, 
very often this also means temporal avaraging [18j |: We 
repeat many times the following two steps: (i) Generate 
the quantum state and (ii) apply a random rotation. The 
disadvantage of this approach is that the execution time 
is proportional to the number of systems in the ensemble. 
Moreover, in many physical realizations of quantum com- 
puting, e.g., in a Nuclear Magnetic Resonance (NMR) 
quantum computer, this approach cannot be used since 
the systems cannot be individually accessed. From the 
numerical point of view, the problem is that averaging 
over the randomly rotated matrices is a very inefficient 
way for calculating the integral in Eq. ((T|). 

Another approach is using group theory to replace the 



2 



integral ([T]) with a sum over a finite number of rotated 
density matrices [19| . This works for small systems and, 
for example, for N = 2 and d = 2 we need to employ 12 
such matrices 0,(13]. However, the number of unitaries 
needed increases rapidly with N and d, making the imple- 
mentation of twirling for large systems difficult this way 
[l6| . Clearly, this approach does not seem to fit ensemble 
quantum computing. From the point of view of a realiza- 
tion on a digital computer, there is the added complexity 
of computing the required unitaries when compared to 
averaging over random unitaries. 

Numerically, there is another approach for replacing 
the integration with a discrete sum. The idea is that 
twirling transforms any quantum state into a U® N in- 
variant state 0] • The density matrix of such a state can 
be written as p — J2k(-^ k )p-^ k where- the Rk basis oper- 
ators are obtained from orthogonalizing the N\ permu- 
tation operators [1, Since twirling does not change 
the expectation values of Rk, we can use these values 
to reconstruct the final Werner state, Pp, on a classical 
computer. However, once more, while this approach is 
feasible for small systems (For N = 2 and TV = 3 qubits 
there are 2 and, respectively, 5 such matrices [1, Q), for 
large ./V the number of permutation matrices increases 
dramatically. 

In this paper we show that a multi-qudit state can be 
approximately twirled by iterating the single non-unitary 
operation 



o+i 



1 r 

Pk+1 - ~ [Pk 



(4) 



where pk are density matrices and Uk are random uni- 
taries. Using this building block, the error of our ap- 
proximation decays exponentially with the number of it- 
erations. The exponent of this decay depends neither on 
the number of qudits nor on their dimension. In con- 
trast, when approximate twirling is realized by averag- 
ing density matrices obtained from multi-lateral random 
rotations, the convergence is algebraic. We show evi- 
dence that the unitary matrices applied can come from 
a highly imperfect source and also demonstrate through 
examples that the random unitaries can be replaced by 
a set of cyclically alternating unitaries, while preserving 
the exponential convergence. Finally, based on the pre- 
vious ideas, we present a quantum circuit for twirling by 
means of controlled unitary gates (see Fig. [JJ . 

Experimental implementation of this operation looks 
feasible in many physical systems. It is important to 
stress that, when applied to an ensemble of many sys- 
tems, our method does not need individual access to the 
individual systems. We show that for a given set of cycli- 
cally alternating unitaries it is possible to obtain general 
statements for the convergence which are valid for all den- 
sity matrices. This makes it possible to design algorithms 
tailored for the operators available in a given physical 
system. Thus twirling can be one of the quantum al- 
gorithms which are especially fitting for realization on a 
quantum computer. On the other hand, when realizing 



- U. 



-a m 



FIG. 1: Twirling can be realized by the repeated application 
of this basic building block where Uk is a random unitary gen- 
erated for the fcth iteration or a unitary chosen from a cycli- 
cally alternating set of unitaries. 'M' represents measurement 
in the computational basis. 



our algorithm on a classical computer, the programming 
and computational effort is extremely small. 

Our paper is organized as follows. In Sec. [Til we discuss 
the usual way twirling is computed on a quantum or a 
classical computer. In Sec. IIIII we present our proposal, 
together also with an analysis of the convergence of the 
approach. In Sec. IIVI we show that our method is quite 
robust against the imperfections of the random number 
generator. In Sec. [V] we show that instead of random 
unitaries, cyclically alternating operators can also be ef- 
ficiently used for twirling. In Sec. IVII we discuss the case 
of large dimensions. In Sec. lVIll we explain how to use our 
ideas for experiments. In Sec. IVIlTI we show how to gen- 
eralize our method for the numerical integration of useful 
formulas over the unitary group. Finally, in Sec. IIXI we 
discuss connections of our research to existing work. 



II. STRAIGHTFORWARD NUMERICAL 
INTEGRATION 

Pp can be approximated by an average of a finite num- 
ber of randomly rotated density matrices 



M-l 



(5) 



fc=i 



Here M denotes the number of terms and we assume that 
the unitaries {Uk} are distributed uniformly in U(d) ac- 
cording to the Haar measure. In this section we examine 
how well PmP converges to Pp for increasing M. 

Since we use random matrices, we will obtain a dif- 
ferent state for each realization of PmP- To analyze the 



3 



error, we introduce an expectation value or average over 



the different choices for Ui. as 21 



(A) 



AdU 1 dU 2 dU 3 



(6) 



Using this average, we can analyze how fast PmP con- 
verges to Pp for increasing number of unitaries. Simple 
calculations show that the average error of a particular 
initial state, p, decreases algebraically as M _1 



(IIPmp-PpII 2 ) 
= (IIPa/pII 2 ) 
= (IIPa/pII 2 ) 
i 



M 



(IMI 



IIPHI 2 
IIPpII 2 

IPpII 2 ). 



2Tr((P M p)Pp) 



(7) 



where ||A|| 2 := Tr(A* A) is the Hilbert-Schmidt norm. In 
the derivation we used that (Pmp) = 0-/M)p + [(M — 
l)/M]PpandTr( j oP j o) = Tr[(Pp) 2 ]. 

While computing the error for a given state is illumi- 
nating, it is more useful to characterize the convergence 
of Pa/ in a manner that is independent of the initial 
state. For that, first we will show how to define a ma- 
trix describing the action of a linear superoperator and 
will define a measure of distance between superoperators. 
Then, we will determine the matrices describing the ac- 
tion of P and P m i and will compute the norm of their 
difference. 

Density matrices are vectors in a Hilbert space of com- 
plex matrices with the scalar product (p, p') := Tr(pp'). 
Thus it is convenient to switch from matrix notation 



P 



= £>w|fc><l| 



hi 



and treat the matrices as vectors defined by [2z 



J2pki\i) ® \k). 



(8) 



(9) 



hi 



That is, p is obtained from p by joining its columns con- 
secutively into a column vector. We can use the vector 
form for any Hermitian operator A, not only for density 
matrices. Then the expectation value of A can be written 
as 



Tr(Ap) = {Afp. 



(10) 



Any physically allowed transformation of the density ma- 
trix is a linear positive map and it can be written as a 
matrix acting on p 



p' = Sp. 



(11) 



Matrix S describes the transformation realized by the 
superoperator. Both the vectors p, p', and the matrix 
S have to satisfy constraints to ensure the Hermiticity 
and the positivity of the density matrices. The distance 



between superoperators can be measured in the form of 
Hilbert-Schmidt norm of their difference |23| 



||S-£|| a := Tr[(S - S)(S - S)t], 



(12) 



In this formalism, the superoperators describing the 
action of P and Pm are, respectively, 



St 



S 



PAI 




M-l 



®2N 



5>* 

fe=l 



(13a) 
,(13b) 



where Id denotes a.dxd unit matrix. Based on Eq. (fl3k ). 
it is easy to sec that 

SpSp = S P = S P . (14) 

Using these, straightforward calculation shows that 



\Spm ^ Spll 2 ) 



= (IIS: 



PAI | 



|S P || 2 -2Tr(Sp(Sp M )) 



= ^(l|lf iV || 2 -||5p|| 2 ). (15) 

Thus the error in the superoperator formalism decays 
algebraically with increasing number of steps, M, irre- 
spective of the initial state (see Appendix |A1 for details). 



III. TWIRLING USING A RECURSIVE 
FORMULA 

In order to decrease the error of the result, rather than 
doubling the number of terms in the summation and com- 
puting Tt*2M, we can apply twice the averaging operation 
with (M— 1) unitaries and calculate PmPmP- In this sec- 
tion we show that, even though in both cases ~ 2M ran- 
dom unitaries are needed, the error of the second method 
is much smaller. 

Let us write out the result after two twirlings explicitly: 



Af-lAf-l 

k=l 1=1 

+ r/ 0Ar o\(u m i f 

^ u M-l+lPlV u M-l+U 

+ {u M -i +l u k r N P [{u M -i + iu k f N ]\ 



(16) 

where {Uk}^ 1 ^ 1 and {Uk}^^ 2 are the random unitaries 
chosen for the first and second twirling, respectively. 
Eq. (TTO|) is the average of M 2 — 1 rotated density matri- 
ces and the original matrix. We have the same number 
of terms when computing P^p. However, the M 2 uni- 
taries are not independent thus we might expect that the 
error for P^/PmP is larger than that for P M 2p. 

Let us now consider repeated applications of P m . For 
simplicity we will first focus on the m = 2 case, leaving 



4 



the m > 2 case for later. After M iterations, the outcome 
is 



Q M p:=P 2 P2...P2/0 = 




Using the definition Eq. (|17p we can write the recursive 
formula 

QmP = \ [Qm-iP + U% N (QM-irftO 1 ] , (18) 

where again Um is a random unitary. As before, we mea- 
sure the convergence of this operator by the average error 
in the Hilbert-Schmidt norm 

(HQmP-PHI 2 ) = (HQmpIP) - ||Pp|| 2 . (19) 

For computing the error as a function of M, we need the 
M -dependence of the (||Qmp|| 2 ) term on the right hand 
side of Eq. (TT9"]) . For that first we express (||QmpI| 2 ) with 

<IIQa/-ipII 2 ) 

(IIQmpII 2 ) = |[<||Q^-ip|| 2 > 

+ (||Qm-ip^ w (Q M -ip)(u® N n 2 )] 

= i((||Q A /-ip|| 2 > + ||Pp|| 2 ). (20) 

Then, from Eq. (|2"0"j) the M-dependence of (||Qmp|| 2 ) can 
be obtained as 

(IIQmpII 2 ) = IIpII 2 + (||Pp|| 2 - ||p 2 ||) (l - 2- M ). (21) 

Substituting Eq. ([21]) into Eq. ((19J) we obtain 

(||Q MP -Pp|| 2 > = (||p|| 2 -||PH| 2 )2- J,/ . (22) 

That is, the squared error decays exponentially with M, 
while according to Eq. ([7]) the decay was oc M _1 for the 
method described in Sec. [TTJ Note that computing QmP 
and PmP need the generation of M and M — 1 random 
unitaries, respectively. Thus the computational effort is 
roughly the same for the two cases. 

One can repeat this calculation in the superoperator 
picture, the The definition of Sqm based on Eq. (fT5|) is 

S QM = 2 { S Q(M-i) + [(Um N )* <8> Ufj N ] Sqfju-i)} • 

(23) 

For the error of the approximation we obtain 

(\\S QM ~ S P \\ 2 ) = (||S qm || 2 ) + ||Sp|| 2 

- (Tr(4 M 5 P )) - (Tr(S QM S P )) 

= (II^qa/|| 2 >-||5p|| 2 . (24) 

For obtaining the M-dependence of the error, we need 
the M-dependence of ||<Sqm || 2 - This is obtained in two 
steps. First, we use Eq. (|23|) to find a recursive relation 
for <||S QM || 2 ) 

(l|5QM|| 2 >-i((||5 Q (M-l)H 2 > + ll^p|| 2 ), (25) 



and then we obtain (\\Sqm || 2 ) without using recursion as 

(IISqmII 2 } = ||lf"|| a + (||5 P || 2 - \\lf N f) (1 - 2-^). 

(26) 

Combining Eq. (|26|) and ((24)) we obtain 

(\\S QM S P \\ 2 ) = (\\lf N \\ 2 \\S P \\ 2 ) 2- M , (27) 

thus the error of the superoperator decays exponentially. 

Let us now consider combining twirl operations 
with more than two unitaries. With that aim we define 

Qk.m := (Pk) M - (28) 

On the one hand, when computing the average error for 
this operator we obtain formulas similar to Eq. ([22|l and 
Eq. (f2~7|) . but with a faster decay oc K~ M vs. the orig- 
inal oc 2~ M . On the other hand, the number of random 
unitaries required increases, as it is now K — 1 per iter- 
ation step. Based on these we can write the dynamics of 
the mean squared error as a function of the number of 
unitaries Njj as 

(WQkm - P|| 2 ) oc exp ^-^Lnu^) . (29) 

Hence one can see that for a given number of unitaries, 
the smallest error is achieved for K = 2. For many ex- 
periments, this is also a good reasoning since the experi- 
mental effort is very often measured in iV^. Thus we will 
consider the K = 2 case in the rest of the paper. 

We have verified numerically the previous results [24j . 
For this we focused on the three-qubit case, for which we 
can compute the twirl operation exactly using the tech- 
niques mentioned in the introduction [See Appendix |A"] . 
In Fig. [2] we plot the error (||SpAf — Sp|| 2 ) averaged over 
10000 trajectories. As the figure shows, the simulation 
results perfectly fit the exponential decay of the error 
calculated theoretically (|2T|. 

IV. SENSITIVITY TO THE IMPERFECTIONS 
OF RANDOM NUMBER GENERATION 

In this section we examine what happens if our random 
number generator does not work perfectly and the ran- 
dom unitaries are not uniformly distributed over U(d). 
An imperfect random unitary generator can be character- 
ized by the distribution f(U) describing the probability 
density for getting U. We will show that if inf u f(U) > 
then our algorithm still converges to the twirled state and 
the error decays exponentially. 

Let us use a simple model for our faulty distribution in 
which with probability p g the unitary is drawn according 
to the probability distribution g(U) while with proba- 
bility (1 — p g ) it is drawn according to the uniform dis- 
tribution. The corresponding distribution function f(U) 
is 

f(U):= Pg g(U) + (l-p g ), (30) 



5 



For obtaining the upper bound in Eq. (|32l) we used 




Tr 



Qm-ipU^ (Q M -i P )(U® N y < ||Qm-ip| 



dWg{W)W® N (Qm-i P )(W^) 



< iiqm-ip|| 2 , 



10 20 30 40 

Number of iterations 



50 



where U is a unitary matrix. From Eq 
dependence of (||Qa//°|| 2 ) / can be deducec 

<ll<Wll 2 }/<IIHI 2 

2 



(33) 
the M- 



\Vp\V-\\p z 



i - 



i+pI 



-M' 



FIG. 2: Mean squared error for the recursive method with 
random matrices when applied on three-qubit states. We plot 
(dotted) the average over 10000 realizations and (solid line) 
the theoretical prediction computed from Eq. (|27[) . For better 
visibility, the error is shown only for every second iteration. 



(34) 



Substituting Eq. (J34J) into Eq. (J3TJ) we obtain 

(hqmp-ppii 2 )/<(iipii 2 -iiphi 2 ) (t^i 

\ 1 ' Pa 



-M 



(35) 

where we used that J UeU(d) dU = L Expectation values The mean square error of the superoperator correspond- 
ing to Qm also decays as oc [2/(l+p g ')]~ AI . Thus we have 
convergence if p g < 1, i.e., if the uniform distribution has 
a non-zero weight in Eq. (f3TJ)) . For functions f(U) that 
satisfy 



>ueu(d) 

over this probability distribution are computed as 

(A) f := I Af(U 1 )f(U 2 )f(U 3 )...dU 1 dU 2 dU 3 ... 
Jueu(d) 

Let us now examine how the usual method described in 
Sec.|n]is affected by such an error of the random number 
generator. We will define by Pm the equivalent of Eq. §5§ 
with our biased probability distribution. The mean value 
of the density matrix obtained from such twirling is 



Mf(U) >0, 



(36) 



+p g (M-l) J dWg{W)W® N p{W® N y 



Here J dW is an integral over the unitary group U(d). 
Taking the limit M — > oo one obtains 



it is always possible to find a decomposition of the type 
Eq. ([3H| such that p g = 1— inff/ f(U). Thus for such prob- 
ability distribution functions our algorithm converges 
and the error decays exponentially. 

Finally, the sufficient condition for the convergence of 
our method Eq. ([55]) can also be formulated for the case 
when f(U) is of the form 



f(U) = f r (U) + J2ck5(U~V k ), 



(37) 



(P MP ) f -» (l- Pg )Pp + Pg j dWg(W)W®" p(W® N )l 

Thus the expectation value of the operator does not con- 
verge to Pp. 

On the contrary, when P2 is applied M times, the 
state of the system still converges to the twirled state 
and the error decays exponentially with M. To show this, 
let us denote the operation above by Qm- As before, we 
measure the convergence of this operator by the average 
error in the Hilbert-Schmidt norm 

(||Q M/ ,-PHr>/-(||QMp|| 2 )/-||Pp|| 2 . (31) 
Hence straightforward algebra yields 



|QmpII 2 )/<^(IIQm-ipII 2 }/ + 



IPHI 



(32) 



where f r : U(d) t—> M, > are constants, 5 is the 
Dirac delta function and are unitaries. In this case 
the algorithm converges if 

inf / f(U)e{U)dU > 0, (38) 

e(U) J(7G(7(<i) 

where for the function e(U) we require that e(U) > and 



V. DETERMINISTIC TWIRLING WITH FEW 
UNITARIES 

The example shown in Sec. If VI demonstrated that even 
if the random unitaries used in our algorithm come from a 
very imperfect source, the algorithm may still converge. 
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on the number of iterations. We look for the c for which 
the decay of the error is the fastest. Through numerical 
optimisation we find that the error of the superoperator is 
the smallest after 50 iterations for c = 1.0894. Fig. El^a) 
shows the results of our numerical calculations for two 
qubits with this value for c. The dashed line shows the 
square of the error for the recursive method using ran- 
dom unitaries described in Sec. IIIIl Note that the error 
for the deterministic method, denoted by disks, decays 
faster than for the random method. 

Let us see a three-qubit example with three cyclically 
alternating unitaries 



Us 

Uy 



A2-K /3(T X 

e i 

z27r/5(7y 

e ) 

__ z27r/3er 2 



(40) 

Fig. [3Jb) shows the results of our numerical calculations. 
Now the error decays somewhat more slowly than in the 
case of the random method. 

The advantage of our approach is that by studying 
the superoperator, we can make general statements, in- 
dependent from the initial state, about the algorithm. 
Thus without a thorough group-theoretical study we can 
show that the error decays exponentially with M and 
the recursive algorithm with the given alternating uni- 
taries can be used for the efficient twirling of two or three 
qubits, respectively. Similar calculations can be carried 
out for several qubits trying out other unitaries or higher 
dimensions can also be investigated. These calculations 
can always consider the gates easily available in an ex- 
perimental implementation. 



FIG. 3: Time dependence of the error (a) for two qubits and 
(b) for three qubits for the deterministic method using two 
and three unitaries, respectively. For better visibility, the 
error is shown only for every second iteration. Dashed line 
indicates the error for the method using random matrices, 
given in Eq. $27} . 

In this section we will examine what happens if these 
unitaries are not random but they are chosen determin- 
istically from a small set such that they are cyclically 
alternating. 

Let us consider the two-qubit case. Common sense 
tells us that we need at least two unitaries since, two 
unitaries are able to generate the elements of U(2) we 
need for twirling [25[. We will see that two unitaries are 
sufficient. 

Let us choose the two unitaries as 

U x := e lCCT *, 

U z := e lc % (39) 

where a x / z are Pauli spin matrices and c is a constant. 
Now we can use the method described at the end of 
Sec. IIIIl to compute the dependence of the superoperator 



VI. TWIRLING FOR QUDITS WITH LARGE 
DIMENSIONS 

In the literature there was a considerable effort to find 
a method for two-qudit twirling which can be realized 
with relatively few quantum gates even if the dimension 
of qudits is large. In this section we show through ex- 
amples that the number of quantum gates necessary for 
two-qudit twirling with our algorithm seems to scale bet- 
ter with the dimension of the qudits than for the algo- 
rithm generating a d-dimensional random unitary uni- 
formly distributed according to the Haar measure. 

A method for generating a random unitary of large 
dimensions was presented in Ref. (26[. The system was 
considered to be a multi-qubit system which fits well for 
many physical realizations. The algorithm presented has 
two steps: (i) Single-qubit random unitaries act on the in- 
dividual qubits. (ii) A nearest-neighbor Ising interaction 
acts on the one-dimensional array of qubits. These two 
steps must be repeated several times. It was found that 
the number of gates necessary for generating a random 
unitary this way scales exponentially with the number of 
qubits n. 

In Ref. it was shown that two-qudit twirling over 
U(2 n ) gives the same result as two-qudit twirling over 
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the Clifford group. This makes efficient twirling possible 
since the number of gates needed for generating a random 
Clifford group element scales polynomially with n. The 
results of Ref. Q were extended to unitary 2-designs in 
Ref. [lj]. 

Let us now examine whether it is possible to find an 
efficient way to realize our algorithm for d > 2. We also 
consider the d = 2™ case. We look for a simple way 
for generating an imperfect random unitary such that it 
still can be used for two-qudit twirling. In particular, we 
would like that the error does not decay slower than when 
using unitarics which are uniformly distributed according 
to the Haar measure. 

We use a slight modification of the algorithm presented 
in Ref. [26j . Our random n-qubit unitary is generated by 
applying first different random unitaries for each qubit, 
then making the system evolve under an fsing Hamilto- 
nian with nearest-neighbor interaction realizing 



exp 



(41) 



where a is a constant and we consider a periodic bound- 
ary condition. Thus we use a single iteration of the 
method presented in Ref. [26(]. The gate requirements 
increase linearly with n for such an algorithm. 

Next we show simulations with the density matrix 
rather than simulations with the superoperator. The rea- 
son is that the size of the superoperator is f 6™ x 1 6™ which 
would make it possible to consider only small systems. 
We calculate the dynamics obtained from our algorithm 
for several random density matrices which have a uni- 
form distribution according to the Hilbert-Schmidt mea- 
sure [27l | . In order to compare trajectories corresponding 
to different density matrices, we compute the normalized 
error 



(llpll 2 -IIPHI 2 ) 



(42) 



It follows from Eq. ([22| that for the method using ran- 
dom matrices uniformly distributed over U(d) we have 

C 1 — 9 — Af 

"norm 6 

Fig. [4] shows the results of our calculations for d = 2 3 
and d = 2 4 . We used a = 1.10 and 1.03, respectively. 
We find that the error decays almost as in the case of 
using random unitaries uniformly distributed over U (d) 
and the difference between the two error curves seems to 
be sub-exponential. 



VII. EXPERIMENTAL REALIZATION 

There are two different situations from the point of 
view of experimental realizations. In many experiments 
several copies of a quantum system are available at a 
time. Very often the systems cannot be individually ac- 
cessed. However, we would like that these systems un- 
dergo different multilateral random unitary rotations. In 
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FIG. 4: Time dependence of the normalized error of the den- 
sity matrix for bipartite twirling for (a) d = 8 and (b) d = 16. 
In both cases 25 realizations are shown. See text for the algo- 
rithm used for generating d-dimensional unitaries. For better 
visibility, the error is shown only for every second iteration. 
Dashed line indicates the error for the method using random 
matrices uniformly distributed over U(d). 



this case, according to our algorithm we have to achieve 
that at each iteration step half of the systems undergo 
a unitary rotation U® , while the other half does not. 
While numerically the mixing of the state p and the ro- 
tated one, U® N p(Uf N )^ is a matter of adding two ma- 
trices, experimentally this mixing can be done with the 
help of a controlled operation. The corresponding quan- 
tum circuit is shown in Fig. [T] For a single step of the 
algorithm we realize P2 . The inputs are the state and an 
ancilla in a superposition state, (|0) + \l))/y/2. The same 
unitary, Uk is applied on all qudits of the state, but only 
when the ancilla qubit is in state 1. Finally, we mea- 
sure the ancilla and the outcome is ±[p + U® k p{Uf Ar ) t ]. 
This basic block is applied M times, each time using ei- 
ther a different, random unitary or a unitary from the 
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finite set as in Sec.fV] Note that the control qubit can be 
a qubit which have a short coherence time compared to 
the other qubits. It rapidly decays to state |0) or |1), and 
can be used as a sort of classical control for permitting 
the unitary rotations on the other qudits. 

In other experiments only a single copy of the state is 
produced. For having an ensemble average of some quan- 
tity, the experiment must be repeated many times. In 
this case our method can be used the following way. Each 
time after the single copy of the state is created, with 50% 
probability we apply U® , then with 50% probability we 
apply Uf N , etc. The number of gates needed in average 
is the half of the number of iterations. If we use the deter- 
ministic version of our method described in Sec. |V] then 
it makes it possible to twirl with a few single-qubit gates. 
This is an advantage in some systems. For example, when 
using photons created with parametric down-conversion 
and post-selection, the single-qubit gates can be realized 
with wave plates. 



VIII. NUMERICAL INTEGRATION OVER U{d) 

As we have already mentioned, our method can be 
used for integrating numerically expressions of the type 
Eq. ([3]). In this section we discuss how to generalize our 
approach for integrating expressions of the type 

I := / Tr{AiU)Tr(A 2 U)...TT(A m U) 

Jueu(d) 

Ti'(BiJ7 t )Tr(B 2 t/ t )...Tr(B„[/ t )d[/, (43) 

where Ak and Bk are d x d matrices. 

Based on the main ideas of the paper, Eq. (|4"3"|) can be 
computed in two steps, (i) First we need to obtain 

M := t U® m <g> {U r f n dU. (44) 
Jueu(d) 

This can be done by iterating the formula 

M k+1 = \[tf m+n) + U® m ® (Ulr n }M k , (45) 

where Mq = 1 and Uk are random unitaries. The series 
Mk will converge very fast to M. (ii) The second step in 
computing Eq. (|43[) is 

I = Tr{M A x ® A 2 ® A m ® Bx® B 2 ® B n ). (46) 

Note that M does not depend on A/, and Bk- Thus when 
we compute Eq. (|4*3"j) for several {Ak} and {Bk}, we have 
to compute M only once. 

These ideas seem to work also when integrating over 
a subgroup of f (d), in particular, over the special uni- 
tary group SU(d). Such integrals appear, for example, in 
quantum chromodvnamics[28l |29| . 



IX. DISCUSSION 

First let us discuss the importance of the fact that our 
algorithm does not require an individual access to the 
systems of the ensemble. This characteristic is important 
since we are presenting the realization of a superopera- 
tor mapping a density matrix to another density matrix. 
Ideally, we want that this mapping works even if the den- 
sity matrix describes an ensemble of very many systems. 
A method which requires an individual access to the sys- 
tems of the ensemble cannot handle this situation. When 
realizing a superoperator in a physical system, it is also 
advantageous that if a pure state is mapped to a mixed 
one then this mixed state is realized as the reduced state 
of a pure state of a larger system [30( . The usual method 
is not able to create such a purification of the output 
density matrix. In contrast, our method can handle a 
very large ensemble. Also, when we apply the quantum 
circuit proposed in this paper for twirling, and we omit 
the measurements then we get a pure state. The twirled 
state is the reduced state of this pure state. 

The algorithm presented in this paper is intimately re- 
lated to other works on random matrices. For instance, 
Ref. (3~H studies the statistical properties of unitary ma- 
trices composed as the product of random unitaries. In 
Ref. [HJ it is proved that the product of a series of ran- 
dom unitaries with nonuniform distribution converges 
exponentially fast to the uniform distribution in many 
cases. 

The relation of our paper to these works is the fol- 
lowing. We also used composed ensembles of unitary 
matrices. However, when looking at Eq. (|16p . we can 
see that our composed unitaries are not independent and 
they are composed from a small set of random matri- 
ces. Thus, especially in the first part of the paper, the 
main goal was to realize twirling on an ensemble of many 
systems using only a few random unitaries, rather than 
realizing twirling using unitaries from an imperfect ran- 
dom source. Note that we assumed that our unitaries 
were drawn from a perfect random source providing uni- 
taries distributed uniformly in U (d). In the second half of 
our paper we found that our results can also be applied 
to the case of an imperfect source or for a deterministic 
algorithm. 

The key point of our algorithm is mixing of a 
subensemble in the original state and the other 
subensemble in which all the systems undergo the same 
multilateral rotation. This mixing can be realized effi- 
ciently both in a classical computer and in a quantum 
computer. Let us analyze the role of mixing pointing 
out something seemingly paradoxical. Let us consider re- 
defining P2 as the application of a unitary U which with 
50% probability it is the identity and with 50% proba- 
bility it is uniformly distributed. That would amount to 
applying a unitary with a distribution 

h{U):= l -5(U-t)+ 1 -. (47) 
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However, if we compute the error of a single application 
of this " new" P2 we find 

(iip 2P -p P ll 2 ), = J \\u® N p(u® N y-p P \\ 2 h(u)du 

= J(\\pf-\\Ppf)h(U)dU 

= ||p|| 2 -||Pp|| 2 , (48) 

which unlike Eq. (|22|) . does not give us an exponential 
convergence. Of course, this is because in the algorithm 
described in this paragraph there is only one component 
present: The application of a random unitary. The other 
component, namely mixing of two subensembles, is miss- 
ing. 

Moreover, while Ref. [13] studied the convergence of 
the distribution of the composed unitaries to the uni- 
form distribution, we studied the convergence of a cer- 
tain quantum operation, namely, twirling. In particu- 
lar, we computed the exponent of this convergence and 
found that it does not depend on N or d. When study- 
ing unitaries composed from random ones with a nonuni- 
form distribution, clearly the requirements for the con- 
vergence of the operator built from these unitaries are 
much weaker than the requirements for the convergence 
of the distribution of the unitaries. It is also easier to get 
general statements on the convergence of the operator for 
a quite wide class of faulty random unitary generators. 
Indeed, we proved that convergence is reached even in 
the case of a very poor random source. 

Note that recursive algorithms can also be applied to 
summing over discrete groups. For example, it has al- 
ready been discussed in Ref. [|| that a random element 
of the Clifford group can be generated by a sequence 
of 0(n 8 ) operations. At each step, with 1/2 probability 
nothing happens, and with 1/2 probability a random ele- 
ment of the generating set is executed. Another example 
is discussed in Refs. (33,|13|. An JV-qubit state can be de- 
polarized by summing over the stabilizer _group [35l l3r| 
of a Greenberger-Horne-Zeilingcr (GHZ, [33]) state or a 
graph state [H, [39[ as 

P = J2 S kPoSl (49) 

k=l 

where {£&} are the group elements. Exploiting that the 
stabilizer group is commutative and that S% — 1, the 
operation Eq. (|39"|) can be realized in N steps. At step 
k with 1/2 probability nothing happens, and with 1/2 
probability gk is executed. Here gk are the N generators 
of the stabilizer group. 

Finally, twirling a completely positive map [161 ] , rather 
than a quantum state, is also a useful procedure. 
Twirling a map makes it possible, for example, to es- 
timate the average fidelity of a physical implementation 
of the map Q . It can be done with a slight modification 
of our algorithm in three steps: (i) Applying the circuit 



Fig. Q] several times with unitaries Ui, U%, Um, (h) ap- 
plying the map, and (iii) applying again the circuit Fig.[T] 

with unitaries U\j, U] vI _ 1 , u\. The control qubit for 

rotations Uk and U k must be the same. 

X. CONCLUSIONS 

Summing up, we have presented a very efficient ap- 
proach for the realization of twirling. Although it is based 
on random matrices, it converges very fast, that is, the 
error decays exponentially with the number of steps dur- 
ing the iteration. Together with the simplicity of the 
method, this means that our algorithm requires very lit- 
tle experimental or computational effort, for the imple- 
mentation in either a quantum or a classical computer. 
We have demonstrated the robustness of the algorithm, 
which converges both in the case of an imperfect random 
number generator and if the unitaries are chosen deter- 
ministically from a small set. In the future, it would be 
interesting to extend our approach to operations which 
realize twirling with a subgroup of U(d). In particular, 
we would like to apply the method described in Section 
IVIIII for integrating numerically over the SU (d) group 
and look for possible applications. 
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APPENDIX A: COMPUTING THE 
SUPEROPERATORS CORRESPONDING TO 
TWIRLING 

For evaluating the right hand side of Eq. (fT5|) , we need 
to know (||Sp|| 2 ). For studying the convergence through 
simulations of individual trajectories we need even the 
matrix Sp. In this appendix, we will study the general 
properties of Sp and determine it explicitly for small sys- 
tems. 

Based on Eq. (Q3|) we can write 

(||S P || 2 ) = Tr(S P 4) = Tr(Sp). (Al) 

Based on Eq. (fbi|) we also know that Sp is a projector 
matrix with eigenvalues and 1, thus (||Sp|| 2 ) must be an 
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integer. We can get to know more about Sp by recalling 
that twirling produces a Werner state and such a state is 
a a linear combination of permutation operators. Since 
these permutation matrices are not linearly independent, 
one has to first orthogonalize them. Let us assume that 
{Rk}k=i are the matrices obtained this way, satisfying 
Ti(RkRi) = Ski where S is the Kronecker symbol. Now 
it is easy to see that we can write the twirled matrix as 



k 

Hence using Eq. (flU)) . we obtain 

S F = ^i? fe (i? fe ) t - 



Thus 



\Spf) =Tr(Sp) = Nj. 



(A2) 



(A3) 



(A4) 



Now, let us determine Sp explicitly for two and three 
qudits. For N — 2 we have Nr — 2 and a possible choice 
of the basis matrices is |7j 



R 1 := 
R 2 ■= 



Id + Vv. 



Id g Id ~ V12 
Vd(d - 1) 



(A5) 



Here Via is the permutation matrix exchanging the two 
qudits and d is the dimension of the qudits. Hence Sp 
can be reconstructed based on Eq. (1A3|) . For N = 3 and 
e? = 2 we have Nr = 5, while for d > 2 we have = 6. 
The basis matrices can be found in Refs. (§, M. 
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